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SYSTEM AND METHOD FOR DYNAMIC 
AEROELASTIC CONTROL 

CROSS-REFERENCE TO RELATED 

APPLICATION(S) 5 

The present application derives priority from U.S. provi- 
sional application Nos. 61/798,078 filed 15 Mar. 2013 and 
61/905,698 filed Nov. 18, 2013. 

to 

STATEMENT OF GOVERNMENT INTEREST 

The invention described hereunder was made in the per- 
formance of work under a NASA contract, and is subject to 
the provisions of Public Law #96-517 (35 U.S.C. 202) in 15 
which the Contractor has elected not to retain title. 

BACKGROUND 

a. Field of Invention 20 

The invention relates to structural monitoring and, more 

particularly, to a hardware and software architecture for 
dynamic modal structural monitoring and control that moni- 
tors the dynamics of a structure such as an aircraft in the 
presence of complex loading patterns using a robust modal 25 
filtering architecture, with integrated distributed shape and 
accelerometer sensing. 

b. Background of the Invention 

Monitoring the states of a structure is useful for safety, 
control and performance optimization in a wide variety of 30 
applications. For example, fuel efficiency and safety are of 
paramount importance to the aviation industry. Toward this 
end, lightweight composite material use in aircraft is now a 
common practice. However, increased flexibility leads to an 
aircraft which is more susceptible to dynamic phenomena 35 
such as gusts, flutter, and buffeting, which negatively affect 
safety and performance and can be potentially destructive. 
Consequently, demand has grown in the aviation industry for 
aircraft that are capable of monitoring and managing such 
dynamics. Typical aeroelastic control methods identify opti- 40 
mal control surface settings for minimum drag from indirect 
sensory data such as wing root load balances, engine RPM or 
thrust/fuel, but such parameters do not give an accurate 
reflection of the true structural state of an aircraft. 

Several approaches have been proposed for more reliable 45 
real-time estimation of structural states. However most are 
frequency based methods which utilize conventional sensors 
to monitor a structure’s natural frequency response. These 
methods are suitable for single-mode vibration suppression, 
but have limited effectiveness for the multi-mode vibration 50 
suppression. Flexible aircraft are highly complex structures 
with multiple vibration modes, and active control of such 
large flexible structures is more readily achieved with mul- 
tiple input multiple output (MIMO) control algorithms. 

An estimator which supports a MIMO type control algo- 55 
rithm is the spatial filter or modal filter. Modal filtering is a 
known spatial filtering technique which transforms physical 
coordinates to modal coordinates. The modal coordinates can 
be used to capture the contribution of both static and dynamic 
deformation in the structure. Modal filtering methods are 60 
generally well known in the art and are disclosed in numerous 
publications including Shelley, S. J., Investigation of Discrete 
Modal Filters for Structural Dynamics Applications, Depart- 
ment of Mechanical and Industrial Engineering, University of 
Cincinnati, 1990. Discrete/continuous modal filters and 65 
dynamic state estimators have typically been used for sensing 
modal coordinates in trusses and bridges. 


2 

A pseudo inverse modal filter, developed by Zhang, was 
implemented by Shelley et al. on a five meter truss structure, 
showing promise for real-time monitoring of the structure. 
Shelley, S., Freudinger, L, Allemang, R. J., and Zhang, Q. 
“Implementation of a Modal Filter on a Five Meter Truss 
Structure.” IMAC IX Conference, Florence Italy, Proceeding 
Volume 2. p. 1036-1044. Unfortunately the pseudo-inverse 
modal filter is not fault-tolerant and fails as soon as one sensor 
fails. 

Shelley et al. later proposed an adaptive modal filter tech- 
nique to compensate for small sensor failure or calibration 
drift. Shelley, S. J., Allemang, R. J., Slater, G. L., Shultze, J. 
F., Active Vibration Control Utilizing an Adaptive Modal 
Filter Based Modal Control Method, 1 \ th International Modal 
Analysis Conference, Kissimmee, Fla., Feb. 1 -4, 1993. How- 
ever, the Shelley et al. adaptive modal filter must be pro- 
grammed with the natural frequency, modal damping, and 
modal residue of the structure. While these parameters can be 
initially calculated, they will tend to drift over time. Shelly’s 
algorithm is also based upon least mean squares (LMS), 
which has similarly bad resistance to gross outliers (large 
sensor bias) as least squares. In flight, the strain can be very 
high. If a sensor failure occurs and the sensor reading is near 
zero while the rest are at normal values, the sensor reading 
near zero is characterized as a gross outlier. This could com- 
pletely bias a modal filter estimate using either least mean 
squares or ordinary least squares. 

Nevertheless, if these past deficits can be overcome the 
modal filter has the potential to benefit the aerospace field 
significantly. What is needed is a system and method for 
dynamic modal structural monitoring and control that inte- 
grates a very large-scale array of sensors including distributed 
strain and acceleration sensing for the purposes of robust 
modal filtering (in lieu of the pseudo-inverse modal filter), in 
real time, while being completely tolerant to multiple sensor 
failures and asymmetric sensor error. 

The present invention provides such a system, and in so 
doing offers more reliable estimates of the modal displace- 
ments and velocities for flexible aircraft structural feedback 
control than the current state of the art, especially in the 
presence of uncertainty and multiple sensor failures. 

SUMMARY OF THE INVENTION 

It is, therefore, an object of the present invention to provide 
a robust and fault-tolerant system for dynamic modal struc- 
tural monitoring and control. 

It is another object to provide a method and system for 
improved aeroelastic control suitable for flight safety critical 
aircraft operation. 

It is another object to implement a real time spatial modal 
filtering system in a robust manner that is insensitive to asym- 
metric sensor noise and sensor failures. 

It is yet another object to implement a practical aeroelastic 
control (where modal frequencies and modal shapes change 
with air speed) that employs real time spatial modal filtering 
for active flutter suppression (AFS), aircraft performance 
improvement, and loads alleviation. 

According to the present invention, the above-described 
and other objects are accomplished by providing a hardware 
and software architecture for dynamic modal structural moni- 
toring that uses a robust modal filter to monitor a potentially 
very large-scale array of sensors in real time, and tolerant of 
asymmetric sensor noise and sensor failures, to achieve air- 
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craft performance optimization such as suppressing aircraft 
flutter, minimizing drag and maximizing fuel efficiency. 

BRIEF DESCRIPTION OF THE DRAWINGS 

5 

Additional aspects of the present invention will become 
evident upon reviewing the embodiments described in the 
specification and the claims taken in conjunction with the 
accompanying figures, wherein like numerals designate like 
elements, and wherein: 10 

FIG. 1 is a block diagram of an embodiment of the present 
system architecture for aeroelastic control and aerodynamic 
optimization. 

FIG. 2 illustrates the FOSS system on wing 2. 

FIG. 3 illustrates the complex X-56A aircraft model of 15 
Example 2. 

FIG. 4 is a schematic of the command shaping module 90 
of FIG. 1 applied configured for design forthe X-56A model 
of FIG. 3. 

FIG. 5 is a diagram of the accelerometer sensor 6 and FBG 20 
sensor 5 control surface locations used for the X-56A state- 
space model of FIG. 4. 

FIG. 6 is a schematic diagram of the uncertain plant and the 
required controller for the X-56A model. This construct may 
be used to test the flight controller 50 of FIG. 1. 25 

FIG. 7 is a schematic of the plant model used for the 
controller design of FIG. 1 to improve rejection of control 
input uncertainty. This construct may be used to test the flight 
controller 50 of FIG. 1. 

30 

DETAILED DESCRIPTION OF THE PREFERRED 
EMBODIMENT 

The present invention is a hardware and software architec- 
ture for dynamic modal structural monitoring that integrates 35 
modal filtering, to monitor an array of sensors using a robust 
modal filter, in real time, immune to asymmetric sensor noise 
and sensor failures. The invention provides a practical, real- 
time approach to monitor the complex dynamics of a flexible 
aircraft using a robust modal filtering algorithm based on the 40 
hybridization of a computationally efficient concentration- 
based estimator and M-estimation. 

Overview 

45 

An embodiment of the present system architecture for 
aeroelastic control is shown in FIG. 1. The plant 10 is a 
mathematical model used to represent the aircraft to be con- 
trolled, which in this case is illustrated as a wing model. The 
wing is equipped with a sensor array 2 for measuring shape 50 
(deflection or strain) at multiple locations on the wing struc- 
ture. 

Sensor array 2 preferably comprises a fiber-optic sensor 
(FOSS) array 2 which is used to measure strain or shape. The 
individual sensors 5 in array 2 can be used by way of modal 55 
filtering to determine the modal displacement states of aero- 
servoelastic state-space models. A distributed array of Fiber 
Bragg Grating sensors is well-suited as FOSS array 2, as 
explained in B. Childers et al., “Use Of 3000 Bragg Grating 
Strain Sensors Distributed On Four Eight-Meter Optical 60 
Fibers During Static Load Tests Of A Composite Structure,” 
Proc. of SPIE, 4332, 133-142 (2001). FOSS array 2 includes 
numerous FBG sensors, placed at regular intervals along 
multi-core (three or more channels in a single fiber) optical 
fibers. A multi-core FOSS array 2 is able to capture defonna- 65 
tion (x,y,z) from its un-deflected position in real time with 
good accuracy, without any knowledge whatsoever of the 


4 

structure to which it is adhered to. The strain measurements 
from each fiber core indicate the direction and amount of bend 
in the fiber. The amount of strain in each core is proportional 
to the bend radius and to the position of each core relative to 
the bend. Alternatively, FOSS array 2 may employ a photo- 
grammetric technique to measure shape at all (x,y,z) loca- 
tions. 

The aircraft carries a computerized feedback loop capable 
of executing the various feedback software modules 10, 20, 
30, 40 of the present invention, including plant model 10, 
modal filter assembler software 20, robust modal filter 30, and 
function approximator 60. Suitable computing devices 
include a processor, memory (e.g. RAM), a bus which 
couples the processor and the memory, a mass storage device 
(e.g. a magnetic hard disk or an optical storage disk) coupled 
to the processor and the memory through an I/O controller. It 
is envisioned that at least a 2.8 GHz processor will suffice to 
provide the necessary computation rate of 1,830 Hzm for 
real-time operation though parallel processing may be pre- 
ferred in some cases. The aircraft also carries a flight control- 
ler 50 which also has various control software modules 60, 
70, 80 for developing and executing the various aeroeleastic 
control commands as a result of the feedback software 10, 20, 
30. 40. 

The modal filter assembler software 20 described herein 
computes the components necessary for the robust modal 
filter 30. Assembler software 20 comprises either of two 
variations. If strain measurements exist in Sensor array 2, 
then the strain is computed at each (x,y,z) location of the FBG 
sensors in real time. Assembler software 20 also assembles 
the strain mode matrix corresponding to the same (x,y,z) 
locations. The measurements and strain mode matrix are 
passed to the robust modal filter 30. Note that the strain mode 
matrix must only be computed once. 

The second variation of software 20 is provided when 
shape at all (x,y,z) locations is estimated from a measurement 
technique such as photogrammetry in real time. During pro- 
cessing of Assembler software 20 the shape modal matrix is 
computed at the same (x,y,z) locations. The measurements 
and shape modal matrix are passed to the robust modal filter 
30. Again, the shape modal matrix must only be computed 
once. 

Getting these distributed measurements and a modal 
matrix (either shape or strain as per above) is all that the 
robust modal filter 30 requires. 

Aircraft aeroelastic behavior is a function of the modal 
coordinates. By controlling modal coordinates, one may con- 
trol the shape of the vehicle. The shape of the vehicle is 
intrinsically related to the aerodynamic performance of the 
vehicle. Thus control of the modal coordinates can lead to 
control of aeroservoelastic performance of the vehicle. The 
robust modal filter software 30 robustly computes the 
required modal coordinates for aeroservoelastic control 
through the computed modal matrices and on-line real time 
distributed sensor (preferably strain) measurements. 

The data from software 20 are input to a robust modal filter 
software module 30, which manages the contribution of the 
thousands of sensors 5 in the system and converts the infor- 
mation through concentration and M-estimation into the 
modal coordinates of the system. Thus, the modal filter 30 
performs a multiple degree of freedom modal coordinate 
transformation. In accordance with the present invention, the 
robust modal filter 30 performs online robust regression 
rather than a least-squares regression, to significantly reduce 
sensitivity to outlier data. Ordinary least squares can give 
misleading results if the error variation of each sensor is not 
normal; thus ordinary least squares is highly susceptible to 
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asymmetric sensor noise and/or sensor failure. The present 
robust regression method (to be described) is designed to be 
not overly affected by violations of assumptions by the under- 
lying data-generating process. A concentration algorithm is 
proposed as a basis for the robust modal filter derivation as 5 
they are applicable to large data sets. E. A. Fox, Q. F. Chen, 
and L. S. Heath. “A Faster Algorithm For Constructing Mini- 
mal Perfect Hash Functions”, Proc. 15th Annual Int’l ACM 
SIGIR Conf. on Research and Development in Information 
Retrieval, pp 266-273 (1992). Typically robust regression 10 
methods for 1 ,000s of measurements are not computable in 
real time. The robust modal filter 30 solves this with a careful 
mix of M-estimation and concentration principles and intro- 
duces a new trim criterion. 15 

The modal coordinates output from the robust modal filter 
30 are input to a Function Approximator 60. Other inputs to 
the Function Approximator 60 may include distributed aero- 
dynamic sensors and optionally fuel gauge information. With 
this information, a functional relationship between the modal 20 
coordinates and the aerodynamic sensors can be generated. 
The relationship is passed to a flight controller 50 which 
includes a diagnostics and performance analyzer module 70. 
The diagnostics and performance analyzer module 70 esti- 
mates performance optimization, stability and control mea- 25 
sures needed to command the shape of the aircraft to an 
optimal shape while staying within constraint boundaries. 
These functional parameters are then sent to a command 
optimizer 80 which confirms the optimal flight system com- 
mands to implement the optimal shape. Finally, the optimal ,(1 
commands are assembled in a command shaping module 90, 
which gives the actual command setpoint to the active control 
system. All the foregoing components are described below in 
more detail. 35 

Sensors 5 

The present invention employs a fiber-optic sensing system 
(FOSS) 2 comprising a plurality of parallel optical fibers 
running along each wing and each having numerous Fiber 
Bragg gratings 5 (“FBGs”) interposed along the length an 40 
optical fiber, spaced at approximately every V 2 ". 

FIG. 2 illustrates the FOSS system 2 on a wing including 
three parallel optical fibers 42, 44, 46 running along wing 2 
and each having a plurality of FBGs 52, 54, 56 . . . interposed 
along the length. To capture sufficient bending information, 45 
the sensors 5 are placed span-wise along the entire wing. To 
capture torsional effects, a plurality of fibers 42, 44, 46 (e.g., 
three) are placed chord-wise. The spacing between each sen- 
sor 52, 54, 56 location along each fiber was set at (A-inch 
intervals. 50 

Light from a tunable laser source (TLS) is split between 
reference and measurement amis of an interferometer. In the 
measurement path, the light is further split to interrogate each 
single-mode optical fiber and return the reflected light. The 
laser is swept so that back reflections from points at different 55 
distances along the fiber under test correspond to different 
beat frequencies. The reflection positions along the fiber can 
be distinguished by their different beat. The software 20 
processes the strain measurements from sensors 5 to derive a 
highly accurate depiction of how each fiber is positioned, and 60 
smooths the data so that derivatives of the modal displace- 
ments can be taken to achieve modal velocities. The assem- 
bler software 20 calculates a modal matrix representing spa- 
tial or strain information at each sensor 5 using the finite 
element method. The assembler software 20 then reduces the 65 
modal matrix according to which sensors 5 are used. The 
distributed measurements from sensors 5 and reduced modal 


matrix (either shape or strain as per above) are then sent to the 
robust modal filter 30 for computation of the modal coordi- 
nates. 

Robust Modal Filter 30 

The design of modal filter 30 from an array of sensors 5 
requires the output signals of each sensor 5 to be weighted 
based on their estimated variance. The measured strain is 
approximately a linear combination of the sensor 5 strain data 
matrix, and so the measured strain is asymmetrically distrib- 
uted. Importantly, wing loading varies with aerodynamic 
condition, and so the underlying strain distribution is difficult 
to predict. Without a known symmetrical distribution, appli- 
cation of most computationally efficient robust outlier detec- 
tion methods is challenging. This strain mode filter problem is 
herein solved using the derived concentrated modal estimator 
(CME) for robust modal filter 30. 

The CME modal filter 30 operates by a recursive sequence 
of steps to provide modal coordinate estimates from data 
which is closest to the computed (in the algorithm step) sta- 
tistical multivariate center of the data. Hie derivation of the 
robust modal filter is comprised of an M-step and concentra- 
tion step. The M-steps occur within each concentration step. 
The CME modal filter 30 concentration step trims outliers, 
and then within this concentration step there are multiple 
M-estimation steps. As described below, during M-estima- 
tion the weights are computed based on the computed vari- 
ance of the sensors 5. The weights used in M-estimation 
downweight or upweight a particular sensor’s importance 
when estimating the modal displacements. Importantly, the 
M-step requires the concentration step (trimming of gross 
outliers based on trim criterion) or it will fall into a local 
minimum and become biased. Therefore the M-estimator and 
concentration estimator of CME modal filter 30 work syner- 
gistically. These substeps are described in more detail below 
and the M-step derivation is given first: 

a. M-Step Derivation 

The CME modal filter 30 utilizes M-estimates within each 
concentration step. M-estimators are well-known gradient 
descent algorithms. They are computationally efficient, affine 
equivariant, robust to masking effects. The strain at measure- 
ment locations may be expanded as a summation of an infinite 
number of orthogonal strain mode shapes (Eq. (0.33)). 


sk(x c , y c , Zc,t) = J^ ?; y c , Zc) 


(7.1) 


To reduce model complexity, only a subset m of mode shapes 
which dominate the response are included in the strain modal 
matrix. It is assumed that the subset of modes captures the 
contributing dynamics and the sensors 5 are subject to ran- 
dom normal error. This introduces an error term into Eq. 
(0.33) which can be modeled as a normal distribution e^eN 
At any discrete time t=t, the quasi-static approximate 
reading of any sensor 5 can be given as in Eq. (0.1), 


Sk(x c , y c , Zc,T) = 2 qi(T)l/li(X c , y c , Zc) + Bk 


(0.1J 



US 9,073,623 B1 


7 


8 


where m is the number of mode shapes retained in the model. 
Consider the linear model for the k r,! sensor measurement to 
be described by Eq. (0.2), 


Sk(Xc,y c ,Zc,T)= (°-2) 

2 q^T)ipi(x c , y c , Zc ) + e k = ¥* (x c , y c , Zc)q(r) + e k 


where e k is a finite residual (i.e. measurement error), ’P J: (x < ., 
y c ,z c )eR lxSr is thek* row of the strain matrix. q(x)eM 4x1 isa 
vector of estimated modal displacements. From the sensor 5, 

6 readings, the obj ective is to estimate q(x) . This can be solved 1 5 
as a maximum likelihood estimation (MLE) problem which is 
posed as minimization of an equally weighted summation of 
a function of the residuals (See Eq. (0.3)). 


b. M-step Operation Within a Concentration Step 
The solution of Eq. (0.5) must be computed within each 
concentration step, c for the proposed concentrated estimator. 
To improve the convergence to the unbiased solution of q(x) 
sensors 5 6 which are most outlying are completely removed. 
For the new group of sensors, M-estimation is used to find 
improved feature estimates. Selection of the influence func- 
tion is critical to the M-estimator’s performance. For its gross 
outlier rejection capability, Tukey’s bisquare (also known as 
biweight) function is the preferred choice. The function is 
chosen to compute the weights with the residuals of the data. 
The bisquare weighting function w=(j)(x)/x is defined for the 
k' ; ' sensor as in Eq. (0.6), 



(0.6) 


Yjp(*k) = J]p(Sk(Xc, y c , Zc, r) -%(x c , y c , Zc)q(r)) 

k=l k=l 

where S is the set of strain sensors and p(x) is an objective 
function with special properties. A reasonable p(x) must be 
symmetric, zero when evaluated at zero, increasing for 
increasing arguments and differentiable. 


where o k <b ’ c> is the median absolute deviation (MAD), h 0 is a 
tuning constant, c is a concentration step, and b is an IRLS 
iteration count of the M-estimator. To achieve the maximum 
95% asymptotic efficiency assuming residuals have a Gaus- 
25 sian distribution, a tuning constant of, for example, h(,=4.685 
is required. 

The MAD for the k ,; ' observation is calculated as in Eq. 
(0.7), 


Define the influence function <(>(x)=p'(x) as the differential 30 
of the objective function p(x). The influence function char- 
acterizes the proportional impact of the residuals on the esti- 
mate. The impact of an OLS residual on the estimate is 
directly proportional to the size of the residual, which is why 
OLS is not robust. To find q(x) the summation given in (0.3) 35 
is differentiated by q(x) and is set equal to zero. By complet- 
ing this, the following equality is achieved (See Eq. (0.4)). 


■i, (0.4) 40 

2j f( s k ( x c , y c ,Zc,T)- % (. x c , y c , Zc)q(r»% (x c , y c , z c ) = o 

4=1 


Let w i .(e i .)=c|)(e it )/e i . for any ^(ej.) then the weighted objective 
function can be rewritten as in Eq. (0.5), 


Yj w * ( e * (*« yc, Zc, T) - % (. x c , y c , Zc)q(r))% ( x c , y c ,Zc) = 0 

*=i 50 


o ! .''’'‘'=MED(le t < i - c> -MED(e ti '- c> )l/0.6745 (0.7) 

where the constant scaling 0.6745 is required to achieve a 
37% Gaussian efficient consistent estimator of the standard 
absolute deviation. While relatively low efficiency, the pur- 
pose of using MAD instead of using the true scale is to resist 
outliers. This it achieves remarkably well since the median is 
high breakdown. However, The MAD is developed for sym- 
metric distributions and does not address distribution skew- 
ness. This may be of concern since the explanatory data is 
multivariate skewed. Improvements of the MAD approxima- 
tion for asymmetric long-tailed distributions are available if 
necessary. Given the weights, vt k < ' b,c) the linear system of 
equations is solved for q (6,c) (x) given sensors in subset S g c as 
in Eq. (0.8). 



( 0 . 8 ) 


y c , Zc, T) - 'i’kiXc, y c , z c )§ <i, ' c) (r))'f J (x e , y c , z c ) = 0 


which results in the weighted least squares problem. Equation 
(0.5) can be solved as a system of equations. Under normal 
conditions an efficient estimate of q(x) can be calculated. The 
weights w^e^.) are affine equivariant and modeled as func- 
tions of the residuals e k . The residuals are dependent upon the 
weights. 

Therefore the method of iteratively reweighted least 
squares (IRLS) is used to solve the optimization problem. 
This proceeds by solving for an initial least squares estimate 
q(x) and computing the residuals and weights. Using the 
weighted observations a new feature estimate q(x) is com- 
puted and the residuals and weights are recalculated. The 
feaUires or modal displacements q(x) of the hyperplane 
approximately satisfying for all sensors, Eq. (0.5) appears 
within a few iterations. 


The weighted least squares problem for the c rh concentration 
step is solved in the same way as Eq. (0.5). Equations (0.6)- 
55 (0.8) are the primary feature estimator equations used within 
the concentration operator. They are iterated within any con- 
centration step for a specified number of M-steps, ^resulting 
in q (6/ ’ c, (x). The next section shows how this estimate is used 
to trim data through concentration before reapplying M-esti- 
60 mation. 

c. Concentration Operation 

The purpose of concentration is to iteratively remove poor 
observations and asymptotically approach the true statistical 
center of the data distribution. Utilizing sensors 5 nearest to 
65 this centroid are assumed to give the best feature estimates. A 
best estimate of this center is the multivariate location T and 
dispersion V of the data. A redescending M-estimator is a 
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robust estimator of multivariate location and dispersion for 
theoretical asymmetric distributions. 

Sensors furthest from this centroid are downweighted in 
Eq. (0.8). However, it is proposed to not just down-weight 
sensors, but to completely remove some gross outliers from 
consideration. The redescending M-estimator does in fact 
equate weights to 0 for gross outliers. The weighted sensor 
removal methodology to improve gross outlier rejection is 
developed here. Let the k sensor data vector be defined as in 
Eq. (0.9), 


10 

tion of D 2 (x a .), k=l ... S is equal to or greater than the 
distribution of D 2 ('P /t (x c ,y £ .,z £ ,)), k==l ... S if the sensor data 
has a Gaussian error distribution. With this assumption, the 
impact of an additional feature may be assumed to change the 
5 distribution of D 2 by the additional degree of freedom (DOF) 
impact in a chi -square distribution. Recall that D 2 is given in 
units of variance. This implies that the variance will increase 
with the additional DOF. Therefore it can be assumed that 
D 2 ( l Fi(x t .,y £ .,z £ .))+0(n J )>D 2 (xj : ), k==l . . . S. Assuming the 
111 adjustment of 0(n^) is due to the noise of the strain data the 
scalar upper bound is defined as in Eq. (0.13): 


x k = PfyUc. yc, Zc)*k(xc, yc , Zc> X)] 


(0.9) 


Defining the data vector in this way ensures that outliers 
which depend on sensor readings shall not occur in the 
explanatory data. Only fixed outliers can occur in the explana- 
tory data, which are to be addressed with a fixed trim crite- 
rion. From any sample sensor set S„°zdS a location vector 20 
(See Eq. (0.10)) 


T (b f ,c) _ 


1 




Vi=l 


(0.10) 

25 


and dispersion matrix [See Eq. (0.1 1)] 


30 


l A b f c ) = 


(b f ,c) T 

V x k x t 


( 0 . 11 ) 


are estimated in the c th concentration step. The weightings sire 
the result of ^iterative M-steps over the subset of sensors S g L 
For the estimated location and variance, the squared 
Mahalanobis distance (D 2 ) (a descriptive statistic that pro- 
vides a relative measure of a data point’s residual distance 40 
from a common point.) is computed for every sensor data 
point k as in Eq. (0.12). 

( 0 . 12 ) 

This multivariate distance differs only from the Euclidean 45 
distance only in that it accounts for correlations between data 
points and is scale-invariant. If the population has a multi- 
variate normal distribution, the D 2 is asymptotically approxi- 
mated by a chi-square distribution. With this knowledge, 
statistical cutoff points from the inverse cumulative distribu- 50 
tion can be determined. The initial distribution of D 2 may be 
computed from the fixed modal matrix and time varying set of 
strain data with Gaussian noise. The maximum of the com- 
puted D 2 may be used as an upper bound for removing gross 
outliers using iterative concentration. During each concentra- 55 
tion step gross outliers are removed and the location and 
dispersion are re-estimated. The sample location and disper- 
sion more closely resemble the population location and dis- 
persion. Therefore the small outliers become more pro- 
nounced. As the D 2 (xj.) increases the sensor can be identified 60 
as an outlier and removed. Outliers missed by this trim pro- 
cedure will more likely be down-weighted in the M-estimate 
[See Eq. (0.8)]. Finding the upper bound D ub 2 is time con- 
suming to implement requiring 1,000s of simulations since 
the strain is time varying. Since most of the data is described 65 
by the constant strain data matrix, an approximation can be 
used for the upper bound. It can be assumed that the distribu- 


te = P c mzx.D 2 (V k (x c , yc, Zc)) 

keS 


(0.13) 


where P c is a tuning constant chosen to be slightly greater than 
1. The tuning constant accounts for 0(n 5 ). By removing a 
portion keS sensors with T> 2 <Y> ub 2 a new candidate group of 
sensors S^ +1 is found for the next concentration step, and 
consecutive M-steps. Simulation studies verify this approxi- 
mation of the upper bound, Y) ub 2 to be good for the strain 
mode matrix and strain data. 

d. Robust Starts and Operations 

Robustness for multi-stage estimators tends to come from 
good starts (initial feature estimates). An initial feature esti- 
mate from a high breakdown estimator is used to start the 
M-estimator for MM-estimates. The median ball algorithm 
proposed by Olive (2004) uses feature estimates from sensors 
closest to the median as a robust start, and is a good choice 
inasmuch as it works reasonably well for nominally skewed 
distributions. Tlie first estimate of the system when x is 0, (i.e. 
when the sensor system is first operational), is calculated with 
a non-robust least squares estimate. Since the first estimate is 
assumed to come from a working sensor system it is a robust 
estimate. The initial robust feature estimate q (o,o:) (0) is found 
by solving the least squares problem presented in Eq. (0.14), 


4 (0.14) 

2 (sk(x c , y c , Zc, T) - 4 ' k (x c , y c , Zc)q ,O ’°'(0))%(x c , y c , z c ) = 0 

k-l 


where S g ° is the set all of the available working sensors. 

During operation, a robust start is paramount. A significant 
advantage of a time based sensor system is that previous close 
estimates are available. The most robust start will therefore be 
the estimate from the previous time step. This is because the 
strain change is assumed to be small between discrete time 
steps. Thus, the robust starts between discrete time steps are 
implemented as in Eq. (0.15), 

^ 0 -°)(t)=^ <4 ^(t-1) (0.15) 

where b^-is the total number of M-steps chosen and c^is the 
total number of concentration steps. 

The importance of starts carries over into the concentration 
steps themselves. In order to be hi gh breakdown, each con- 
centration step requires a robust start. Since the initial start 
q(°’ o, (0) is robust the final estimates at the end of each of the 
concentration steps: 

q (A/,1) (x), q (6/,2) (x), . . . q (A/ ’ c/_1; (x) are robust. This mean that 
the estimates of corresponding concentration steps are robust 
starts for respective next concentration steps: 
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q(°’ 2, (V), q (0 ’ 3) (t), . . . q (0 ’ c/) (T). Therefore the following inlier- 
itance rule [See Eq. (0.16)] is used to generate robust starts 
between concentration steps: 

x)=f b f-‘\x) (0.16) 

The full steps of the CME modal filter 30 for any discrete 
time step are summarized in Algorithm 1 assuming that an 
initial OLS feature estimate has already been computed with 
Eq. (0.14) at time 0. 

{s k (x c ,y„z c ,%)4 b W(x-l)}^ b f- c i > Algorithm 1 

1) If c==0, compute q (0,0) (T) with Eq. (0. 1 5). Else compute 
q (0 ’ c \x) with Eq. (0.16). 

2) For b=0: by, Iteratively compute weights, vi ] } bf ’ d) with 
Eqs. (0.6)-(0.8) with q (0>c) 

3) Compute location T (See Eq. (0.10)) and dispersion V 
(See Eq. (0. 1 1 )) with computed w s (4/, ' 3) 

4) Compute D 2 (x t ), k=l . . . S, using Eq. (0.12) with 
computed T and S 

5) Generate a new sensor set S^" 1 by trimming sensors 
below cutoff D ub 2 described by Eq. (0.13). 

6) If c<Cy, Go to Step 1. Else output c^ bfiCf \x). 

For each time step, the M-step iteration count by may be 
initialized to be large, so that a robust redescending M-esti- 
mate initializes the CME modal filter 30. This improves the 
algorithm’s stability during the concentration steps. After- 
wards, single M-steps where byis equal to 1, may be utilized. 
This has the effect of improving computational efficiency. 

The CME modal filter 30 given by Algorithm 1 is robust 
estimator with more than 0.5 breakdown point given the use 
of robust starts which are estimates from previous time steps. 
By assuming that the initial estimate at time 0 occurs when the 
sensor system is operating normally, the OLS estimate at time 
0 is robust. This robust estimate is used as a robust start to the 
next time step; the corresponding robust estimate is used for 
the start in the next time step and so on. Therefore robustness 
is guaranteed between time steps. Robustness is also guaran- 
teed between concentration steps. Therefore, iterative appli- 
cation of the robustness inheritance concept guarantees 
robustness through all concentration operations. 

Function Generator 60 

The CME modal filter 30 outputs the modal coordinate 
estimates as an input to the Function Generator 60. The air- 
craft must also provide the network with local aerodynamic 
information and possibly fuel gauge information. The Func- 
tion Generator 60 takes the modal coordinates from CME 
modal filter 30 and aerodynamic information and iteratively 
updates a functional relationship between the two. 

The function generator 60 approximates an unknown func- 
tion of the modal coordinates F(modal coordinates) based on 
aerodynamic data, such that the aerodynamic state=F(modal 
coordinates). Of course, in order to determine what F is, local 
aerodynamic data must be provided and so some measure- 
ment of aerodynamic state is necessary. Consequently, aero- 
dynamic state sensor(s) such as hot film strips or leading edge 
sensors for monitoring pressures over the wing surface (for 
optimal drag configuration), fuel rate sensors (for fuel effi- 
ciency), or any other suitable sensor(s) for measuring aero- 
dynamic state are provided. For example, given fuel rate 
sensors it is possible to form the functional relationship fuel 
rate=F (modal coordinates). Once the function generator 60 
approximates the function F, it becomes possible to know 
where to command the modal coordinates to in order to 
achieve the optimal aerodynamic state (e.g., optimal drag 
efficient shape). 

The function generator 60 is passed to the flight controller 
50 for analysis. From the analysis, a minimum drag aircraft 
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control surface configuration for fuel savings can be attained, 
or flutter suppression, etc. In effect, aircraft inflight shape 
optimization becomes possible given the structural states 
measured. The function approximator 60 may comprise a 
5 neural network that employs on-line learning and computes 
functional relationship changes in the plant 10. The Flight 
Controller 50 computes the required control surface deflec- 
tions for an optimal drag configuration. An example of a 
aeroelastic control neural network is presented in U.S. Pat. 
10 No. 7,617,166 to Haudrich et al. 

Flight Controller 50 

The function generator 60 is input to a Flight Controller 50 
in the vehicle, which may be an existing flight controller 
programmed with the requisite control modules 70, 80, 90 of 
15 the present invention. Specifically, Flight Controller 50 
includes a diagnostics and performance analyzer module 70 
which analyzes the function output from function generator 
60 and estimates the optimal modal coordinates to correct to, 
as necessary to achieve performance optimization, stability & 
20 control measures needed to command the shape of the aircraft 
to an optimal shape. The optimal modal coordinates from 
performance analyzer module 70 are sent to an Optimizer 80, 
which is a function approximator similar to 60 and can be 
interrogated to find the optimal control surface settings. Opti- 
25 mizer 80 may likewise be a neural network or a less-complex 
function solver such as a genetic algorithm. A genetic algo- 
rithm is a known method for finding globally optimal solu- 
tions fora given function. See, e.g., Mitchell, “An Introduc- 
tion to Genetic Algorithms”, Cambridge, Mass., MIT Press 
30 (1996). Finally, the optimal control surface settings from 
Optimizer 80 are sent to Command Shaping Module 90 
which computes the optimal flight system commands to 
implement the optimal shape and compiles the command set 
for input to the existing active control system. 

35 The following examples illustrate the above-described sys- 

tem in more detail. 

Example 1 

40 The above-described system has been verified with a 
simple FEM model of a wing developed using MATLAB/ 
Simulink™. The wing model allowed the user to input 
aspects such as wing span and chord, as well as number and 
location of control surfaces. The finite element model was 
45 coupled with an unsteady aerodynamics model and gust 
model captured by Rational Function Approximations 
(RFAs) of the generalized forces found through frequency 
based linear aerodynamic methods. Fighter aircraft type 
actuators were utilized for dynamic control surface simula- 
50 tion. State space models were formed and a simulated modal 
filtering system was designed. The state space models were 
also utilized for control design. The system was implemented 
with an array of fiber optics with FBG sensors. The wing 
model was structurally simple. 

55 

Example 2 

The above-described system has been verified with a com- 
plex X-56A aircraft model as shown in FIG. 3. The design of 
60 the Command Shaping Module 90 (from FIG. 1) is shown in 
FIG. 4. This architecture assumes a virtual deformation con- 
trol architecture with modal reference tracking. Since defor- 
mation is directly related to modal displacements, the signifi- 
cant modal displacements are commanded directly, thereby 
65 leading to deformation changes in the vehicle. This is what is 
referred to by virtual deformation control architecture. Note 
that the FOS simulation and FOS Noice Simulations at right 



US 9,073,623 B1 


13 

would be replaced by the FOSS Sensing System 2 of FIG. 1 in 
a real-aircraft implementation. 

Identification of significant modes in an aircraft structure is 
similar to that for a wing model. V-g and V-f plots or equiva- 
lent may be used to determine interacting modal coordinates, 5 
and modal coordinates must be selected for feedback which 
contribute to flutter, modal vibration and are within the actua- 
tor bandwidth. Selection of these modes can be an iterative 
process. Once feedback modes have been decided upon, the 
state space matrices are updated. The matrix equation of to 
motion of the aeroelastic system in discrete coordinates is 
given as in Eq. (0.17). 

[M hl J{q}+[C hh ]{q}+[K hl J{q}+[M h J$}= 

(0.17) 

Analogous to Euler’s first law, F=ma (or in this case ma=F), 
the left side of the aeroelastic system equation represents the 
structure’ s mass properties and applied motion while the right 
side represents the aerodynamics forces. The aero module 
includes several components such as aero paneling, aero 2 o 
modes and gust modes, as well as GAF calculation and RFA 
design. Points between the FE model and aero model are 
selected that accurately represent the force transferal from the 
aerodynamic control points to the structural grid points. This 
allows calculation of the GAFs from modal deformation at 2 s 
structural grid points. It is necessary to convert the frequency- 
domain GAF to the Laplace domain. However, the Laplace- 
dotnain unsteady aerodynamics must be in a rational function 
form to be incorporated into the time-domain state-space 
equations of the aeroelastic system. Therefore, the unsteady 30 
aerodynamic forces are approximated through the following 
RFA as in Eq. (0.18). 
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ZAERO™ software system by Zona Technologies, Inc., a 
software tool for advanced aeroelastic design and analysis (by 
running open-loop ASE analysis in frequency domain first, 
e.g., normal flutter analysis, and then comparing with open- 
loop analysis in the time domain using state-space models to 
ensure the RFA is accurate). The state vector is shown in Eq. 
( 0 . 21 ), 


*(r) = x g {t), **(0, x hg (t), jwWI 7 (0 ' O1) 

where x AF (t)6R^' s<1 is a vector of airframe states, x e (t)e 
R , " ><1 is a vector of modal displacements, x 5 (t)eR mxl is a 
vector of modal velocities, x^tleR^ 1 is a vector of aero- 
dynamic lag states and x acf (t)e]R is formed from a vector 
of actuator accelerations, velocities, and displacements. 

Trim conditions were applied to the nonlinear equations of 
motion to develop the state space models. The trim analysis 
ensures gravity is included into the model to provide more 
accurate flight responses for rigid body motion. In addition, a 
1G gravitational force was applied and the only free variable 
was the angle of attack. These conditions lead to a trim angle 
of attack and a trimmed elastic solution. Therefore all the state 
space models’ states are considered to be perturbations about 
this point. 

The state space models must also model actuator dynam- 
ics. A third order transfer function of the actuators is formu- 
lated for each control surface to incorporate the actuator 
dynamics into the system. The actuator transfer function is 
represented in Eq. (0.22). 


[fit!*,)] * [Aoi + [/4il(£fe;) + [A 2 ](iki) 2 + 


(0.18) 35 


<hP) an 

u aci CO -s 2 + a a s 2 + a a s + a i3 


( 0 . 22 ) 


[D][»/[/i- j[«i] ‘rat/*,-) 


From principles of analytic continuation, the unsteady aero 40 
modeled in terms of reduced frequencies can be expanded 
into the entire Laplace domain through the following substi- 
tution given in Eq. (0.19), 


The LTI transformation resolves the transfer function into 
actuator state space matrices of the as in Eq. (0.23). 


Hac,)-- 


o 1 

0 0 


0 1 tot 

1 {Xac,} + \ 0 U 


(0.23) 


—at 3 ~ a ,2 -at i 


V 


s = g + ik 


(0.19) 


45 


where s is the Laplace Variable and g is the non-dimensional 
damping. Aero lags are used to model the delayed force 50 
effects of the unsteady aerodynamics. The aero lags are mod- 
eled as zeros along the negative real axis of the Laplace 
domain since the aero force translation to the structure is not 
real time. By setting g=0, for the non-dimensional damping, 
the aerodynamic forces are now expressed in the Laplace 55 
domain as in Eq. (0.20). 


[ew] = 


Qhc\ 


= 1AJ 


( 0 . 20 ) 


■p[tiih + ■pjl'tih 2 + 


[D][d/l-^[«l] Kb 


These models are appended to the state space models for all 
ten control surfaces, increasing the state space order by the 
order of the actuator transfer function multiplied by ten. 

Sensors 5 which measure strain were also defined from 
modal matrix computations in NASTRAN™ (a finite element 
analysis (FEA) program originally developed for NASA in 
the late 1 960s). The shape modal matrices were converted to 
strain modal matrices. The known simulation modal coordi- 
nates are used to back out strain measurements assuming 
superposition for small strain. This was done in lieu of using 
actual fiber optic sensor strain measurements, which would 
replace the above methodology. 

The modification of the state space matrices for modal 
coordinate feedback is very similar to how the wing model 
state space matrices are modified. Consider the relationship 
of the strain-based modal filter to the modal displacement 
state vector x^(t) as in Eq. (0.24), 


Roger’s method was chosen to compute the RFA. Con- 65 
straints are applied at specific reduced frequencies (such as at 
zero) and aerodynamic derivatives may be matched using the 


x s (t)=q(t)=( x V T ’ < ¥- lx V r s rn (x c ,y c ,z c t)+£ (0.24) 

where ’P={S li (|> I (x < _,y £ ,,z c ),i=l ... N} is a NxN strain matrix 
with N strain modes, s m (x c ,y c ,z c )€l A><1 is the measured 
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strain and 6 is a normal error. This represents the modal filter 
in its pseudo-inverse form. Since the modal filter gives a 
partial state estimate of the full state vector, x(t) as defined in 
Eq. (0.2 1 ) the form of the output matrix is simply identity for 
measured modes. The sensor output matrix, C, is a matrix of 5 
row vectors relating the sensory information to the state vec- 
tor, x(t). The output matrix is adjusted for rigid body state 
feedback concurrent with modal deformation state feedback. 

It is assumed that m modes are retained for measurement. The 
output matrix is formed assuming all modal displacements, to 
x s , and airframe states, x^p, are measured, as shown in Eq. 
(0.25). 


IaSxAS 

OmxAS 


OASxm 0/1 Sxm 0 AS xf 
tmxm 0 m y m 0 m y J 


0/lSxz 
0 myz 


(0.25) 15 


As with the wing model, there is no requirement to measure 
all flexible states (or rigid states) to adequately sense the 20 
vehicle state. Higher frequency modes generally do not sig- 
nificantly contribute strongly to the deformation and may be 
ignored to reduce the size of the C matrix. 

After selecting and modeling modal coordinates, the state 
space matrices can be reduced in order and used to design the 25 
control laws. Similar to the wing model (Example 1 ), control 
design can proceed with any desired control methodology. 

For the wing model, the p-Optimal design methodology was 
chosen. See, S. Skogestad and I. Postlethwaite, Multivariable 
Feedback Control Analysis and Design, West Sussex, 30 
England: John Wiley & Sons, Inc, 2005. Once the controller 
has been designed which meets requirements, the modal 
coordinates are also fixed. They can be used in the modal filter 
design. 

The modal filter 30 for the X-56A is in units of strain 35 
measured directly from the FOSS sensors 5. To utilize a strain 
mode matrix for modal filter 30, the axial strain modes must 
be computed. The close proximity of the sensors (of approxi- 
mately one-half-inch) along the SFOS gives a unique oppor- 
tunity to calculate axial strain modes directly from the defor- 40 
mation at the sensor locations. An approximation of the strain 
mode matrix was made using the computed deformations for 
each mode shape at the sensor 5 locations. This relies upon a 
simple strain to displacement relationship which is a deriva- 
tive. This can also be done by fitting a shape polynomial to the 45 
surface and taking a derivative at the desired points. 

The transformed mode shapes defined at sensor locations 
are collected to form the modal matrix. These mode shapes 
are pure elastic and thus are used to form a proper modal 
reference transformation. The modal reference command is 50 
shown in Eq. (0.26): 

q n j(t)=<&*(l n m r )d n .j(x c ,y c ,z c ,t) (0.26) 

where d re ^(x c ,y c ,z c ,t) is a vector of deformations from the 
un-deformed aircraft corresponding to the index vector, \ r , in 55 
the pure elastic deformation modal matrix, O. 

Recall that each row of the modal matrix corresponds to a 
physical location on the aircraft, (x c ,y c ,z c ). The index of 
modes, m,„ corresponds to the index of modes within the 
modal matrix, which the controller is designed to track. The 60 
selection of 1,. is a research topic within itself. However, here 
it is recommended to select points with high deformation in 
the commanded mode shapes. A better technique would also 
select points which minimized residual mode amplitudes. 
This completes the modal filter 30 design methodology. The 65 
following sections put the presented modal filter design meth- 
odology into practice on the X-56A vehicle. 


Once strain estimator 20 has captured the modal displace- 
ments of sensors 5, the Modal Filter 30 computes the modal 
matrix or transformation matrix, which is expressed as a sum 
of modal vectors. 

The X56A is intentionally designed to encounter three 
flutter instabilities within its operating flight envelope: body 
freedom flutter (BFF), symmetric wing bending torsion flut- 
ter (SWBT), and antisymmetric wing bending torsion flutter 
(AWBT). BFF is a phenomenon where the rigid body mode of 
pitch and plunge couples with the elastic mode of first wing 
bending. The other two modes are traditional elastic flutter 
modes. 

The BFF mode is demonstrated by the interaction of the 
scaled states, rigid body heave state, bz cm , and modal dis- 
placement 8SWB1. The flutter characteristics of the BFF 
mode are divergent oscillations of contributing modes at the 
same frequency, out of phase. 

The SWBT mode also illustrates the unfavorable coupling 
of modal displacements: 8SW1B and 8SW1T. The nature of 
the coupling is difficult to discern because it appears random 
in nature. This means that there is an interaction with another 
mode, likely the BFF mode. 

The AWBT mode shows an in-phase interaction of modal 
displacements: 8AW1B and 8AW1T. The scaled modal 
amplitudes are small but clearly grow in time, verifying that 
this flutter mode also exists at the design flight condition. 
Without control, the aircraft will experience strong flutter 
instability and will require AFS at this speed. 

The flutter modes indicate that four flexible modal dis- 
placements should be used for control feedback, including 
8SW1B, 8SW1T, 8 AW IB and 8 AW IT. However, the fre- 
quencies of the anti-symmetric modes are prohibitively high 
and close to the actuator bandwidth. Therefore only the sym- 
metric E' bending and 1" torsional modal displacements are 
selected for feedback. The requirements of the vehicle con- 
troller are now overviewed, and a p-optimal controller is 
designed, which can track both rigid body states and the first 
two modal displacements. 

Since flutter is a potentially destructive phenomenon, the 
modal controller must be robust to uncertainty. To reduce 
these uncertainties, the following uncertainty requirements 
were defined: 10 percent multiplicative uncertainty on the 
inputs and outputs, and 1 0 percent additive uncertainty on the 
scaled plant. In addition, the controller must be robust to a 
3-percent speed variation. The controller must also achieve 
robust performance despite the aforementioned uncertainty 
conditions. The uncertain plant 10 and required control sys- 
tem 50 may be tested using the construct summarized in FIG. 
6. This represents the requirements which the controller 50 
must meet, and not necessarily how it is designed. Note that 
for other aircraft, these requirements may change. 

The robust stability condition which must be met for the 
p-optimal controller is given as in Eq. (0.27). 

RS <=> (i(M A (;(o))<l,Vo) (0.27) 

The (j. is calculated over the frequency range with the relation 
shown in Eq. (0.28), 


l4.M\( ito)) = 


1 

min!/:,,, | det(/ - k. m M X \U,,)!S\ = 0 
for structured A, cr(A) < 1} 


(0.28) 


where k m is the stability margin, a is the maximum singular 
value, and A is the structured uncertainty. The transfer func- 
tion matrix from the input of the uncertainty blocks to the 
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outputs of them as shown in Figure Error! Reference source 
not found, is presented as in Eq. (0.29), 


M a 


-W,T, -W,KSo -W,KSo 


W*S, 


-W a KS 0 


-W a KS 0 


L GeGv 


W 0 GS, W 0 S 0 -W 0 T 0 


5 

(0.29) 


10 


where W 7 is a matrix of proper input weights, W A is a matrix 
of proper additive weights, W G is a matrix of proper output 
weights, and K is the controller. From this matrix, the salient 
sensitivity and complementary sensitivities from the M A 
structure are identified which correspond to the requirements. 
The magnitude of the singular values predict the performance 
of the control system. These closed loop transfer functions are 
defined as in Eq. (0.30), 

20 

Tj=KG(I+KG)~ l 

KS 0 =K(I+GK)~ l 

25 

S 0 =(I+GK)~ l 

Sj=(I+KG)~ l 

, 30 

T 0 =GK(I+GK)~ l (0.30) 


where T 7 is the input complementary sensitivity, S 0 is the 
output sensitivity, S 7 is the input sensitivity, and T 0 is the 
output complementary sensitivity. 

To improve rejection of control input uncertainty, the con- 
troller 50 was designed around the plant construct shown in 
FIG. 7. 

The design plant, G des , at the design speed, V des , was of 
order 130. To improve controller synthesis G des was scaled by 411 
the full range of actuator movement and expected sensor 
output changes. The translational states, bx cm , hy cm , and 
8z cm , and the velocity state were removed from the model. 
With the resulting 126th-order model, balanced reduction 
was performed to bring the model order to 90 states. The 
selected states to be tracked were: SW1B and SW1T modal 
displacements, pitch angle, 0, and bank angle, 4>. The angle of 
attack a, angle of sideslip |3, and yaw angle, were also 
sensed but were chosen to be suppressed. 50 

Traditional proper weights from a mixed H„ synthesis 
were utilized along with multiplicative uncertainty at the 
plant inputs. The input uncertainty weight, W 7 , was adjusted 
to achieve maximum amplitude near the actuator break fre- 
quency. Sensitivity weight, W s was adjusted for integral 55 
tracking on tracking states and for suppression on suppres- 
sion states. 

It was found that the break frequencies of the modal dis- 
placement performance weights had to be increased 10 rad/s 60 
relative to the airframe state weight break frequencies of 1 
rad/s. The break frequency of the control weight, W„, was 
adjusted to 5 rad/s for reduced control surface movement. The 
break frequency of the complementary sensitivity weight, 
W r , was set to 30 rad/s to improve high frequency noise 65 
rejection. The uncertainty transfer function, N A , was calcu- 
lated to be as shown in Eq. (0.31): 


W,T, 

-WjKSo 

-W s GS, 

-W s So 

-W U T/ 

-W U KS 0 

W t GS i 

- W t T 0 . 


Thus, by reducing the FK norm of N A , the specified uncer- 
tainties in FIG. 6 can be rejected. 

For this control design architecture, the p-optimal control- 
ler was computed using MATLAB’s Robust Control Tool 
Box. To find the controller, a DK-iteration was performed as 
per S. Skogestad and I. Postlethwaite, Multivariable Feed- 
back Control Analysis and Design , West Sussex, England: 
John Wiley & Sons, Inc. (2005), which solves the iterative 
optimization problem shown in Eq. (0.32). 



(0.32) 


The DK-iteration resulted in a 1 62th-order controller after 
some trial and error with the weights in FIG. 7. The controller 
was then internally balanced and truncated to 44 states with- 
out a substantial loss of robustness or performance. This 
resulted in an H, norm of 3.29. 

A concentrated modal estimator (CME) was then devel- 
oped to form the basis for CME Modal Filter 30. The CME is 
primarily a concentration algorithm with re-descending 
M-estimates used in place of OLS within the concentration 
steps. The CME utilizes a fixed trim criterion and a more 
robust start to address asymmetry of the strain mode matrix 
data. For the X-56A sensor system fourteen modal displace- 
ment features are included in the strain mode matrix tp cor- 
responding to 1 ,530 sensor entries. Therefore, the problem is 
computationally burdensome and also multivariate. The 
strain at measurement locations may be expanded as a sum- 
mation of an infinite number of orthogonal strain mode 
shapes (See Eq. (0.33)). 


Skitc, yc,zc,') = Yj ?i(0M*c. yc,zc) 


(0.33) 


To reduce model complexity, only a subset m of mode shapes 
which dominate the response are included in the strain modal 
matrix. It is assumed that the subset of modes captures the 
contributing dynamics and the sensors are subject to random 
normal error. This introduces an error term into Eq. (0.33) 
which can be modeled as a normal distribution e^eN^i^a,,. At 
any discrete time t=x, the quasi-static approximate reading of 
any sensor can be given as in Eq. (0.1), 


Sk(x c , yc, zc, T )=/, qMifrAxc, yc, zc ) + e k 


( 0 . 34 ) 


where m is the number of mode shapes retained in the model. 
Consider the linear model for the k r/ ' sensor measurement to 
be described by Eq. (0.2), 
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Sk(x c ,yc,zc,r) = (0.35) 

y, ?,(x) i// : (x c , yc , Zc) + e t = % (x c , yc, Zc)q(r) + e k 

i=l 5 

where e k is a finite residual (i.e. measurement error), ’P lt (x c ,, 
y^zJeR lxk is the k* row of the strain matrix. q(x)eR 4x1 is a 
vector of estimated modal displacements. From the sensors 
readings, the objective is to estimate q(x). This can be solved 
as a maximum likelihood estimation (MLE) problem (See 
Huber [195]) which is posed as minimization of an equally 
weighted summation of a function of the residuals (See Eq. 
(0.3)), 15 

s s (0.36) 

2 p(e * ) = 2 pO* ( x c , yc.zc.r)- % (x c , yc, zc)q(T)) 


where S is the set of strain sensors and p(x) is an objective 
function with special properties. A reasonable p(x) must be 
symmetric, zero when evaluated at zero, increasing for 
increasing arguments and differentiable. 25 

Define the influence function 4>(x)=p'(x) as the differential 
of the objective function p(x). The influence function char- 
acterizes the proportional impact of the residuals on the esti- 
mate. To find q(x) the summation given in (0.3) is differenti- 
ated bye q(x) and is set equal to zero. By completing this, the 30 
following equality is achieved (See Eq. (0.4)). 


s (0.37) 

2 (*c, yc. Zc, t) - T* (. x c , y c , Zc)q(T)Wk (xc. yc, Zc) = 0 35 

*=i 


Let w k (e k )=ty(e k )/e k for any 4>(e A .) then the weighted objective 
function can be rewritten as in Eq. (0.5), 

40 


20 

of the data. The bisquare weighting function w=<j)(x)/x is 
defined for the k r/ ' sensor as in Eq. (0.6), 


„Je\ ri Jfl -f- 
1 UH r w 


< ho 


otherwise 


(0.39) 


where o k h ' c> is the median absolute deviation (MAD), h 0 is a 
tuning constant, c is a concentration step, and b is an IRLS 
iteration count of the M-estimator. To achieve the maximum 
95% asymptotic efficiency assuming residuals have a Gaus- 
sian distribution, it lias been shown that a tuning constant of 
h 0 =4.685 is required. 

The MAD for the k r,< observation is calculated as in Eq. 
(0.7), 

Ojt (i..o =M ED(le t (i, -''-MED(e (i '-‘’')l)/0.6745 (0.40) 

where the constant scaling 0.6745 is required to achieve a 
37% Gaussian efficient consistent estimator of the standard 
absolute deviation [249]. While relatively low efficiency, the 
purpose of using MAD instead of using the true scale is to 
resist outliers. This it achieves remarkably well since the 
median is high breakdown. 

However, The MAD is developed for symmetric distribu- 
tions and does not address distribution skewness. This may be 
of concern since the explanatory data is multivariate skewed. 
Improvements of the MAD approximation for asymmetric 
long-tailed distributions are available if necessary. Given the 
weights, w k ib ' c) the linear system of equations is solved for 
q (i,c )(x) given sensors in subset S c as in Eq. (0.8). 


(0.41) 

X l ^' C, (^)( i * (jEc ’ ^ Zc ’ ~ 


%(*c, yc, Zc)q <b ' c \Tjftk(x c , yc, zc) = 0 


s (0.38) 

2j w k ( e k )(Sk (xc, yc, ZC, x) - 
*= 1 

% (Xc. yc. Zc)9( T))% (Xc, yc, Zc) = 0 45 

which results in a weighted least squares problem. Equation 
(0.5) can be solved as a system of equations. Under normal 
conditions an efficient estimate of q(x) can be calculated. The 50 
weights v/ k (e k ) are affine equivariant and modeled as func- 
tions of the residuals e k . The residuals are dependent upon the 
weights. Solve for an initial least squares estimate q(x) and 
compute the residuals and weights. Using the weighted obser- 
vations a new feature estimate q(x) is computed and the 55 
residuals and weights are recalculated. The features or modal 
displacements q(x) of the hyperplane approximately satisfy- 
ing for all sensors, Eq. (0.5) appears within a few iterations. 

The solution of Eq. (0.5) must be computed within each 
concentration step, c for the proposed CME estimator. To 60 
improve the convergence to the unbiased solution of q(x) 
sensors which are most outlying are completely removed. For 
the new group of sensors, M-estimation is used to find 
improved feature estimates. Selection of the influence func- 
tion is critical to the CME estimator’s performance. 65 

For its gross outlier rejection capability, Tukey’s bisquare 
function is chosen to compute the weights with the residuals 


The weighted least squares problem for the c ,h concentration 
step is solved in the same way as Eq. (0.5). Equations (0.6)- 
(0.8) are the primary feature estimator equations used within 
the concentration operator. They are iterated within any con- 
centration step for a specified number of M-steps, by-resulting 
in q (6 /’0(x). 

The purpose of concentration is to iteratively remove poor 
observations and asymptotically approach the true statistical 
center of the data distribution. The foregoing estimate is used 
to trim data through concentration before reapplying M-esti- 
mation. Utilizing sensors nearest to this centroid are assumed 
to give the best feature estimates. A best estimate of this 
center is the multivariate location T and dispersion V of the 
data. Sensors furthest from this centroid are downweighted in 
Eq. (0.8). However, it is proposed to not just down-weight 
sensors, but to completely remove gross outliers from con- 
sideration. The redescending M-estimator does in fact equate 
weights to 0 for gross outliers. Let the k* sensor data vector be 
defined as in Eq. (0.9), 


Xk = 1 ^k(Xc, yc, Zc)s k (x c , yc, Zc, x)l (0 ' 4 ’ 

Defining the data vector in this way ensures that outliers 
which depend on sensor readings shall not occur in the 
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explanatory data. Only fixed outliers can occur in the explana- 
tory data, which are to be addressed with a fixed trim crite- 
rion. From any sample sensor set S„ C =»S a location vector 
(See Eq. (0.10)) 


22 

estimate. Since the first estimate is assumed to come from a 
working sensor system it is a robust estimate. The initial 
robust feature estimate q (O ’°\0) is found by solving the least 
squares problem presented in Eq. (0.14), 







M 

w k J x k 


(0.43) (0.47) 

^ (sk (Xc, yc, Zc, r) - % (x c , yc, Zc)q ( °’°\ 0))% ( x c , yc, Zc) = 0 

*=i 

10 


and dispersion matrix [See Eq. (0.1 1)] 




(0.44) 15 


are estimated in the c‘ h concentration step. The weightings are 
the result of iterative M-steps over the subset of sensors S g c . 20 

For the estimated location and variance, the squared 
Mahalanobis distance (D 2 ) is computed for every sensor data 
point k as in Eq. (0.12). 

D 2 (x k )=x k -I b f^)(V b fi cy )-\x k -T b f' ci ) T (0.45) 25 

This multivariate distance differs only from the Euclidean 
distance only in that it accounts for correlations between data 
points and is scale-invariant. If the population has a multi- 
variate normal distribution, the D 2 is approximated and sta- , f) 
tistical cutoff points can be determined. The initial distribu- 
tion of D 2 may be computed from the fixed modal matrix and 
time varying set of strain data with Gaussian noise. The 
maximum of the computed D 2 may be used as a upper bound 
for removing gross outliers. During each concentration step , 
gross outliers are removed and the location and dispersion are 
re-estimated. The sample location and dispersion more 
closely resemble the population location and dispersion. 
Therefore the small outliers become more pronounced. As the 
D 2 (x ft ) increases the sensor can be identified as an outlier and 4(| 
removed. Outliers missed by this trim procedure will more 
likely be down-weighted in the M-estimate [See Eq. (0.8)]. 

An approximation is used for the upper bound D ui 2 . It can be 
assumed that the distribution ofD~(x /t ), k=l ... S is equal to 
or greater than the distribution of D 2 (*P 4r (x ( ,,y ( ,,z ( ,)), k=l . . . S 
if the sensor data has a Gaussian error distribution. Recall that 
D 2 is given in units of variance. This implies that the variance 
will increase with the additional DOF. Therefore it can be 
assumed that D 2 ('P i .(x £ ,,y £ ,,z c ))+0(rL)>D 2 (xj : ), k=l . . . S. 
Assuming the adjustment of 0^) is due to the noise of the 5Q 
strain data the scalar upper bound is defined as in Eq. (0. 1 3): 


D lb = PfOiax D 2 ^ t (.v c - yc, Zc )) 

keS 


(0.46) 

55 


where is a tuning constant chosen to be slightly greater than 

1. The tuning constant accounts for 0(n,). By removing a 
portion keS sensors with D 2 -D„ 6 2 a new candidate group of 
sensors S^ +1 is found for the next concentration step, and 60 
consecutive M-steps. Simulation studies verify this approxi- 
mation of the upper bound, D ub 2 to be good for the strain 
mode matrix and strain data. 

Robustness for multi-stage estimators tends to come from 
good starts (initial feature estimates). The first estimate of the 65 
system when x is 0, (i.e. when the sensor system is first 
operational), is calculated with a non-robust least squares 


where S^° is the set all of the available working sensors. 

During operation, a robust start is paramount. A significant 
advantage of a time based sensor system is that previous close 
estimates are available. The most robust start will therefore be 
the estimate from the previous time step. This is because the 
strain change is assumed to be small between discrete time 
steps. Thus, the robust starts between discrete time steps are 
implemented as in Eq. (0.15), 

? (0 -°'(t)=^/-54(t-1) (0.48) 

where b^-is the total number of M-steps chosen and c^is the 
total number of concentration steps. 

The importance of starts carries over into the concentration 
steps themselves. In order to be high breakdown, each con- 
centration step requires a robust start. Since the initial start 
q <0,0) (0) is robust the final estimates at the end of each of the 
concentration steps: q (6/,1) (x), q ibfi2 \x), . . . q (A/ ’ C ' rl) (x) are 
robust [205], 

This mean that the estimates of corresponding concentra- 
tion steps are robust starts for respective next concentration 
steps: q (0,2) (x), q (0,3> (x), . . . q (0,c/:) (x). Therefore the following 
inheritance rule [See Eq. (0.16)] is used to generate robust 
starts between concentration steps: 

<7 10 -" 1 Xx)=q (b Ffi(x) (0.48) 

The frill steps of the CME for any discrete time step are 
summarized in Algorithm 1 assuming that an initial OLS 
feature estimate has already been computed with Eq. (0. 1 4) at 
time 0. 

{s k (x c ,y c ,z c ,T)4 (l *- c f>(T - 1 ) } Algorithm 1 

1 ) If c=0, compute q (0,0) (x) with Eq. (0. 1 5). Else compute 
q (0 ’ c) (x) with Eq. (0.16). 

2) For b=0: b y. Iteratively compute weights, w 4 5 A/ ’ c ' ) with 
Eqs. (0.6)-(0.8) with q (0,c) 

3) Compute location T (See Eq. (0.10)) and dispersion V 
(See Eq. (0.1 1)) with computed w jt (A/,<2) 

4) Compute D 2 (x^), k=l . . . S, using Eq. (0.12) with 
computed T and S 

5) Generate a new sensor set S g c_i by trimming sensors 
below cutoff D ub 2 described by Eq. (0T3). 

6 ) If c<c f , Go to Step 1. Else output q fA/,c/) (x). 

For each time step, the M-step iteration count b^ may be 
initialized to be large, so that a robust redescending M-esti- 
mate initializes the CME. This improves the algorithm’s sta- 
bility during the concentration steps. Afterwards, single 
M-steps where byTs equal to 1, may be utilized. This has the 
effect of improving computational efficiency. The robust esti- 
mate from previous concentration steps are passed as robust 
starts to respective next concentration steps. Therefore, itera- 
tive application of the robustness inheritance concept guar- 
antees robustness through all concentration operations. 

It should now be apparent that the above-described system 
integrates modal filtering and distributed sensing, to monitor 
an array of sensors using a robust modal filter, in real time, 
immune to asymmetric sensor noise and sensor failures. The 
invention provides a practical, real-time approach to monitor 
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the complex dynamics of a flexible aircraft using a robust 
modal filtering architecture and a computationally efficient 
concentration-based and M-estimator-based robust algo- 
rithm. 

It should be understood that various changes may be made 
in the form, details, arrangement and proportions of the com- 
ponents. Such changes do not depart from the scope of the 
invention which comprises the matter shown and described 
herein and set forth in the appended claims. 

The invention claimed is: 

1. A computerized system for controlling aeroelastic 
behavior of an aircraft, comprising: 

a fiber-optic sensor array comprising a plurality of optical 
fibers oriented substantially parallel along a wing of said 
aircraft, each said optical fiber having a plurality of 
spaced fiber Bragg Gratings (FBGs); 
a processor in communication with said fiber-optic sensor 
array, said processor including a non-transitory com- 
puter-readable storage device storing a plurality of soft- 
ware modules including, 

a plant model comprising a finite element model of said 
aircraft, 

an assembler software module comprising a series of 
instructions for performing the steps of calculating a 
reduced modal matrix representing one of spatial or 
strain information at each FBG in said fiber-optic 
sensor array, 

a modal matrix software module for computing modal 
coordinates, 

a function estimator for estimating an optimal aerody- 
namic state from said modal coordinates; and 
a flight controller for estimating a performance optimiza- 
tion adjustment to said modal coordinates to maintain 
said optimal aerodynamic state. 

2. The computerized system according to claim 1, wherein 
said modal matrix software module is a concentrated modal 
estimator (CME) configured to recursively estimate said 
modal coordinates from said FBG data that is closest to a 
statistical multivariate center of data. 

3. The computerized system according to claim 1, wherein 
said flight controller is configured to compute required con- 
trol surface deflections to produce an optimal drag configu- 
ration for said aircraft. 

4. A system for dynamic aeroelastic control of an aircraft, 
comprising: 

a sensor array attached to said aircraft, said sensor array 
comprising a plurality of fiber Bragg Grating (FBG) 
sensors; 

a processor in communication with said sensor array, said 
processor including a non-transitory computer-readable 
storage device storing software comprising a series of 
instructions for performing the steps of, 
maintaining a finite element model (FEM) of said aircraft 
comprising a set of finite elements interconnected at a 
plurality of nodes, 

mapping a spatial position of each of said FBG sensors, 
monitoring strain at each of said FBG sensors in real time; 
calculating a sensor strain matrix representing linear axial 
strain at each ofsaidFBG sensors from said modal shape 
deformation; 

removing rows of sensor strain matrix based on which 
sensors are used; 

applying a robust modal filter software module for modal 
coordinate transformation of said reduced sensor strain 
matrix into a strain-based modal matrix by projecting 
said reduced sensor strain matrix from physical space to 
modal space to define projected modal coordinates; and 
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a function estimator for estimating the functional relation- 
ship of modal coordinates and aerodynamic efficiency; a 
controller for performance optimization adjustment 
with said modal coordinates to maintain said optimal 
5 aerodynamic state. 

5. The system according to claim 4, wherein said aircraft 
comprises a wing and said is plurality of fiber Bragg Grating 
(FBG) sensors are spaced along a plurality of optical fibers 
running substantially parallel along said wing, 
to 6. The system according to claim 4, wherein said FEM 
comprises a wing model. 

7. The system according to claim 4, wherein said FEM 
comprises a model of the entire aircraft. 

8. The system according to claim 4, wherein said step of 
15 monitoring deformation at each of said FBG sensors com- 
prises measuring strain at each of said FBG sensors and 
estimating modal deformation from said measured strain. 

9. The system according to claim 4, wherein said step of 
reducing said sensor strain matrix comprises weighting data. 

20 10. The system according to claim 9, wherein said step of 

weighting data includes zero-weighting outlier data. 

1 1 . The system according to claim 4, wherein said control- 
ler comprises a flutter suppression module for minimizing 
aircraft flutter. 

25 12. The system according to claim 4, wherein said control- 

ler comprises a drag minimization module for minimizing 
aircraft drag. 

13. The system according to claim 4, wherein said control- 
ler comprises a performance optimization module for maxi- 

30 mizing aircraft fuel efficiency. 

14. A system for dynamic aeroelastic control of an aircraft, 
comprising: 

a sensor array attached to at least a portion of said aircraft; 

a processor in communication with said sensor array, said 
3 5 processor including a non-transitory computer-readable 
storage device; 

a plurality of software modules each comprising a series of 
instructions stored on the machine-readable storage 
device and embodying executable code run by the pro- 
40 cessor to perform a series of steps, said software mod- 
ules further comprising, 

an assembler module for mapping a spatial position of each 
sensor in said sensor array, for monitoring deformation 
at each said sensor, and for calculating a sensor matrix 
45 representing spatial position and strain modes at each of 
said sensors in said array; 

a robust modal filter software module for modal coordinate 
transformation of said reduced sensor strain matrix into 
a strain-based modal matrix. 

50 15. The system according to claim 14, wherein said soft- 

ware modules further comprise a finite element model (FEM) 
comprising a set of finite elements interconnected at a plural- 
ity of nodes. 

16. The system according to claim 14, wherein saidrobust 
55 modal filter software module accompli shes modal coordinate 

transformation by projecting said sensor strain matrix from 
physical space to modal space to define projected modal 
coordinates. 

17. Hie system according to claim 16, wherein saidrobust 
60 modal filter software module accomplishes modal coordinate 

transformation by the substeps of, 

projecting sensor data from said sensor array onto the 
nodes of said FEM using a robust regression operator to 
define a data matrix, 

65 reducing said data matric, and 

projecting said reduced data matrix from physical space to 
modal space to define projected modal coordinates. 
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18. The system according to claim 14, further comprising a 
function estimator for estimating an optimal aerodynamic 
state from said modal coordinates. 

19. The system according to claim 18, further a control 
module for estimating a performance optimization measure 5 
to be imparted to a control surface to maintain said optimal 
aerodynamic state. 

20 . A method for dynamic aeroelastic control of an aircraft, 
comprising the steps of: 

providing said aircraft with a sensor array; to 

maintaining a finite element model (FEM) of said aircraft 
comprising a set of finite elements interconnected at a 
plurality of nodes, 

mapping a spatial position of each of said FBG sensors, 
monitoring deformation at each of said FBG sensors from 15 
an un-deflected position in real time; 
calculating a sensor strain matrix representing spatial posi- 
tion and strain modes at each of said FBG sensors from 
said deformation; 

reducing said sensor strain matrix; 20 

applying a robust modal filter software module for modal 
coordinate transformation of said reduced sensor strain 
matrix into a strain-based modal matrix by projecting 
said reduced sensor strain matrix from physical space to 
modal space to define projected modal coordinates; 25 
estimating an optimal aerodynamic state of said aircraft 
from said projected modal coordinates; and 
estimating a performance optimization adjustment to said 
modal coordinates to maintain said optimal aerody- 
namic state. 30 



